#ifndef __NURBS_MESH__
#define __NURBS_MESH__

#include "attribute.h"
#include "primitive.h"
#include "node.h"

class Element;

class Mesh : public Primitive {
public:
	virtual Node* GetNode(size_t i) = 0;
	virtual const Node* GetNode(size_t i) const = 0;
	virtual Element* GetElement(size_t i) = 0;
	virtual const Element* GetElement(size_t i) const = 0;
	const std::vector<Node*>& GetBoundary() const {return boundary;}
	const std::vector<Node*>& GetBoundary() {return boundary;}
	
	template<typename T>
	void AddNodeAttribute(const std::string& attribute, const T& t) {
		std::vector<T> x(SizeNodes(), t);
		nodeAttribute.Add(attribute, x);
	}
	
	template<typename T>
	void AddElementAttribute(const std::string& attribute, const T& t) {
		std::vector<T> x(SizeElements(), t);
		elementAttribute.Add(attribute, x);
	}

	const Attribute& GetElementAttribute() const {return elementAttribute;}
	Attribute& GetElementAttribute() {return elementAttribute;}

	const Attribute& GetNodeAttribute() const {return nodeAttribute;}
	Attribute& GetNodeAttribute() {return nodeAttribute;}
	
	virtual size_t SizeNodes() const = 0;
	virtual size_t SizeElements() const = 0;

	bool ExistNodeAttribute(const std::string& attribute) const {
		return nodeAttribute.Exist(attribute);
	}
	
	bool ExistElementAttribute(const std::string& attribute) const {
		return elementAttribute.Exist(attribute);
	}
protected:
	Attribute nodeAttribute, elementAttribute;
	std::vector<Node*> boundary;
};

void Transform(const std::string& attribute, Mesh* mesh, const Mat4x4n& T);
void TransformConstrained(const std::string& attribute, Mesh* mesh, const Mat4x4n& T, const std::string& markAttribute, const size_t mark);

template<class T>
std::vector<Node*> FindNodes(Mesh* mesh, const std::vector<T>& P, const std::string& attribute) {
	std::vector<Node*> node;
	for(size_t i = 0; i < P.size(); i++) {
		nfloat min = std::numeric_limits<nfloat>::max();
		size_t kmin;
		for(size_t k = 0; k < mesh->SizeNodes(); k++) {
			const T& Pk = mesh->GetNode(k)->Get<T>(attribute);
			const nfloat l = length(Pk - P[i]);
			if(l < min) {
				min = l;
				kmin = k;
			}
		}
		node.push_back(mesh->GetNode(kmin));
	}
	return node;
}

template<class T>
void LaplaceSmooth(Mesh* mesh, const nfloat lambda, const std::string& source, const std::string& destination, const std::string& boundary) {
	for(size_t i = 0; i < mesh->SizeNodes(); i++) {
		Node* v = static_cast<Node*>(mesh->GetNode(i));
		if(v->Get<size_t>(boundary) == 0) {
			nfloat lambdaj = lambda;
			const T& pin = v->Get<T>(source);
			T& pout = v->Get<T>(destination);

			T sum = (T)0.0;
			for(size_t j = 0; j < v->NumNeighbors(); j++) {
				sum = sum + (v->Neighbor(j)->Get<T>(source) - pin);
			}
			sum = (1.0 / v->NumNeighbors()) * sum;
			pout = pin + lambdaj * sum;
		} else {
			v->Get<T>(destination) = v->Get<T>(source);
		}
	}
}

template<class T>
void LaplaceSmooth(Mesh* mesh, size_t steps, nfloat lambda, const std::string& destination, const std::string& temp, const std::string& boundary) {
	std::cerr << "Laplace Smooth..." << std::endl;
	for(size_t i = 0; i < steps; i++) {
		std::cout << (int)(100.0 * ((nfloat)i/(steps - 1))) << '\r';
		LaplaceSmooth<T>(mesh, lambda, destination, temp, boundary);
		LaplaceSmooth<T>(mesh, lambda, temp, destination, boundary);
	}
}

template<class T>
void ZeroNodes(Mesh* mesh, const std::string& attribute) {
	for(size_t i = 0; i < mesh->SizeNodes(); i++) {
		Node* n = mesh->GetNode(i);
		n->Get<T>(attribute) = T(0.0);
	}
}

template<class type, size_t dim>
math::Vec<type, dim> Center(const std::string& attribute, const Mesh* mesh) {
	math::Vec<type, dim> center;
	for(size_t i = 0; i < mesh->SizeNodes(); i++) {
		const Node* n = mesh->GetNode(i);
		center += n->Get< math::Vec<type, dim> >(attribute);
	}
	return (1.0 / (mesh->SizeNodes())) * center;
}

template<class type>
void Clamp(Mesh* mesh, type a, type b, const std::string& attribute) {
	type max = -std::numeric_limits<type>::max();
	type min = std::numeric_limits<type>::max();
	for(size_t i = 0; i < mesh->SizeNodes(); i++) {
		const type& s = mesh->GetNode(i)->Get<type>(attribute);
		if(s > max)
			max = s;
		if(s < min)
			min = s;
	}
	
	type scale = 1.0 / (max - min);
	for(size_t i = 0; i < mesh->SizeNodes(); i++) {
		Node* v = mesh->GetNode(i);
		type& s = v->Get<type>(attribute);
		type u = (s - min) * scale;
		s = u * (b - a) + a;
	}
}

#endif //__NURBS_MESH__
